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ABSTRACT 

In this paper we investigate the polarization property of the radiation amplified by 
■ astronomical masers in the presence of a strong magnetic field. Our model explicitly 

takes into account the broadband nature of the radiation field and the interaction 
' of the radiation with the maser transition J = 1—0. The amplification of different 

realisations of the background continuum radition by the maser is directly simulated 



and the Stokes parameters of the radiation field are then obtained by averaging over the 
+3 , ensemble of emerging maser radiation. For isotropic pumping and partially saturated 

masers we find that the maser radiation is linearly polarized in two representative 
cases where the magnetic field B makes an angle 6>=30° and 8—90° to the maser axis. 
The linear polarization for maser radiation obtained in our simulations for both cases 
are in agreement with the results of the standard model. Furthermore, no instability 
■ during amplification is seen in our simulations. Therefore, we conclude that there is 

£SJ , no problem with the previous numerical investigations of maser polarization in the 

C*~) ■ unsaturated and partially saturated regime. 
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1 INTRODUCTION 



k>( \ Strong maser radiation has been detected in various astronomical environments, from nearby starforming regions to distant 
« I ■ galactic nuclei. The extremely high brightness temperature of astronomical masers allows very high angular resolution ob- 
CS ' servations to study the kinematics and physical properties of the masing environments (Reid & Moran 1981). It has been 
recognized that maser radiation is generally polarized, either linearly or circularly. In some cases the polarization degrees 
up to 100% have been observed. Because polarization of the radiation field is related to spatial anisotropy, i.e the presence 
of a strong and ordered magnetic field, maser radiation holds the promise of providing valuable information on the elusive 
magnetic field. 

A complete understanding of the mechanism to generate polarization of maser emission is essential to interpret observational 
results and to infer the physical properties, such as magnetic field direction and strength, in the maser environment. Theoret- 
ical study of astronomical masers started with the work of Litvak (1970). Subsequently Goldreich et al. (1973) worked out the 
basic model of maser polarization for the J=l-0 transition. Their work predicts that, in the presence of a strong magnetic field, 



radiation from saturated masers is linearly polarized, reaching a limit of Q/I=— 1 for sin 9 ^ | and (3sin 9 — 2)/3sin 9 
for sin 2 9 ^ |. However, the result is derived only for the case of a fully saturated maser and the question of how that 
limiting polarization is reached for any real astronomical maser remains open. Numerical studies by Western & Watson 
(1984), Deguchi & Watson (1990) and Nedoluha & Watson (1990) provide the dependence of fractional linear polarization 
on the saturation parameter R/F, the ratio between the stimulated emission rate R and the loss rate V. Numerical results 
show that for the J=l-0 transition, the linear polarization increases slowly toward the limiting solution, and for higher lying 
rotational transitions such as J=2-l, the limiting polarization predicted by Goldreich et al. (1973) can only be attained in 
the fully saturated regime, namely R/F" ~ 100 or even higher. Observationally, such high intensities are almost impossible to 
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be realised in astronomical masers. More detailed calculations by Deguchi & Watson (1990), taking into account population 
mixing between magnetic sub-levels in the regime of large stimulated emission rate, actually predict that the polarization 
fraction decreases and disappears at high enough intensity. Thus, the limiting polarization can never be reached in high lying 
transitions such as J=2-l. 

Observationally, SiO maser lines in vibrationally excited state v = 1, such as J=3-2 and even J=5-4, from evolved stars have 
been known to possess large linear polarization (Mcintosh & Predmore 1993). This discrepancy is often cited (Elitzur 1993, 
1995) as deficiency of the framework under which previous numerical studies are carried out. Recently, in a series of papers 
(Elitzur 1991, 1993, 1996) a different model of the maser polarization was proposed. The linear polarization fraction is deter- 
mined through finding the eigenvalues of the radiation transfer equations. In an important departure from previous studies, 
an ensemble average of the different modes of the radiation field, which is inherently broadband and random, is considered. 
Some novel solutions are found, namely that unsaturated maser emission can have the same polarization fraction as saturated 
masers and when the magnetic field is close to the maser axis (sin 2 6 ^ |), propagation of polarized emission is forbidden. In 
addition, the same analysis predicts that polarization properties of astronomical masers are spin independent, i.e high lying 
transitions should behave in the same way as the often studied J=l-0 transition. Instability during the amplification of maser 
radiation is cited as the main reason for the new solutions, although in subsequent publications (Elitzur 1995, 1996), more 
emphasis is put on the ensemble average over modes or configurations of the radiation field. Elitzur (1995) points out that a 
full understanding of the creation and evolution of polarization in masers requires simulations involving a statistical ensemble 
of waves. These new results are vehemently criticized by Watson (1994), who performs again the stability analysis on the 
radiation transfer equations and finds no instability. More recently, Gray (2003) points that the multi-level model (Field & 
Gray 1988, Gray & Field 1995) can be formally reduced to the idealized two-level case with a similar set of equations as 
in Watson (1994), leading to the same predictions for the linear polarization properties of astronomical masers as obtained 
earlier by Watson and co-workers. 

Obviously, the explicit incorporation of the broadband random radiation field into the study of astronomical masers is of 
geat importance. That kind of simulation might help to elucidate the current debate regarding the polarization theory of 
astronomical masers. To simulate from first principles the amplification of radiation by a masing medium, as suggested by 
Elitzur (1995), is a complicated task. In addition to a new formulation of the interaction process between radiation field and 
the maser medium, a dramatic increase in computing power to follow the evolution of a large number of realisations of the 
incident radiation field is also required. 

So far the only work aimed to treat the broadband maser radiation from first principles, and in a transparent and self- 
consistent manner, is published in Menegozzi & Lamb (1978). However, due to limited computing power, they could not 
perform a large enough number of simulations to effectively draw any firm conclusion on the statistics of the maser radiation 
field. In a previous paper (Dinh-V-Trung 2009) we have followed the formulation of Menegozzi & Lamb (1978) to investigate 
the standard theory of a scalar maser and the statistics of the maser radiation field. In this paper we generalize our model of 
the astronomical maser to include the vector nature of the radiation field. That will allow us perform simulations of polarized 
maser radiation, which is the main focus of this paper. We hope that our work will provide a small step toward a better 
understanding of the properties of astronomical masers. 



2 BASIC THEORY 
2.1 Radiation field 

In this paper we consider a one-dimensional maser and assume that a strong magnetic field exists inside the maser medium 
and makes an angle 9 to the propagation direction of the radiation. That means the Zeeman splitting gcuB between magnetic 
sub-levels of a particular rotational level J is much greater than the decay rates due to collisions and/or radiative transitions. 
As a result, the magnetic field defines a preferred direction in space and a good quantization axis. We adopt here the geometry 
used by Deguchi & Watson (1990) in which the magnetic field B is aligned with the Oz' axis of the B-frame. The k-frame 
which has the Oz axis directed along the propagation direction is obtained from the B-frame system by rotating through an 
angle 6 about the Ox' axis. 

The spherical basis is defined in the conventional way (Brink & Satchler 1994, Zare 1988): 
e p =±i = =F (e x + e y ) /a/2; 

e p =o = e z (1) 
where (e x ,e y ,e z ) are the unit vectors along the frame axis. Similarly the spherical bases in the B-frame is defined as: 
e M =±i = =F (e x / + e y / ) /\/2 
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The relation between the two bases associated with the k-frame and the B-frame follows easily: 
e p= ±i = =F(e x ' ± icos#e y / ± i sin# e z /)/\/2 
e p =o = sin6*e y / + cos#e z / 



(3) 



We represent the electric field of the maser radiation and the induced macroscopic polarization vector in the spherical bases: 



E(z,t) 



Re 



£(-l) P E p e_ p 



P(z,t) = Re 



(-l) P P P e- 



(4) 



Where E p and P p are spherical components of the electric field and polarization vector, respectively. 

Because our aim is to simulate the spectral properties of the maser radiation, it is preferable to work from the beginning in 
the frequency domain. In expressing the frequency dependence of the electric field and polarization vector, we retain only the 
positive frequencies. The anti-resonant (negative) frequencies are ignored. This approximation is usually referred to as the 
rotating wave approximation. The spherical components of the electric field and polarization vector can be written as: 



Ep(z,t) 
PpOM) 



E p (z,t)e-' u ' o(t - z/c) 
P p (z,t)e-' u ' o(t - z/c) 



(5) 



Where luq = 2nvo is the angular frequency at the maser line center. The amplitudes E p (z,t) and P p (z,t) are assumed to 
vary slowly with time in comparison to the term e~ 1 " '. Using the spectral representation theorem (Priestley 1981), for any 
realisation of the radiation field during the time interval T the amplitude of electric field and polarization vector can be 
expressed in terms of Fourier series: 



E p (z,t) = E p {z,uj n )e 



iu> n (t — z/c) 



P p (z,t) 



n— + 00 

E 



P p (z,u„) e 



-iw n (t — z/c) 



(6) 



where u) n = 2im/T. The Stokes parameters of the radiation field are the quantities measured directly in observations of 
astronomical masers. These parameters are commonly used to characterize the polarization properties of the maser radiation. 
Because the radiation field is stationary and ergodic, these parameters can be defined as either an ensemble average or time 
average of different realisations of the radiation field. For the convenience of presenting the simulation results, we will follow 
the same convention as in Deguchi & Watson (1990) in defining the four parameters {X, Q,U, V) for each frequency uj n during 
an interval T of each realisation: 



[Q(Wn) 



T(u) n ) Auj 
V(w n )Aw 
HJ(lo„)] Aui 



[E-{u n )E*_(u n ) + E+(b}„)El(u n )] 
[E-(u n )E*_ (w„) - E+{u n )E* + {u n )] 



8tt 
c 

8^ 



(7) 



In this definition, all the parameters (J, Q,U,V) are real, for example Q and U can be obtained directly from the real and 
imaginary part of the last equation. The inclusion of Alo—2tt/T is due to the fact that we use Fourier series to decompose 
the radiation field into discrete harmonic components. Each Fourier component represents the radiation field in the frequency 
band of Au). The coherence time of the electric field (Mandel & Wolf 1965) is indeed the interval T. Thus, during this interval 
of time, the radiation field within the frequency band Aw can be considered as quasi-monochromatic. The four parameters 
(I, Q,U, V), which are real and defined as above, are similar to Stokes parameters of a monochromatic wave and statisfy the 
usual relation: 



1 (w„) 



g 2 K)+wV„) + v 2 K) 



(8) 



We emphasize here that these parameters contain all the information on the amplitudes and the relative phase between 
different spherical components (E + and E-) of the radiation field in each frequency band. The usual Stokes parameters of 
the random radiation field are then the ensemble averages, denoted as (...), of the above (I, Q,U, V) parameters: 



J(w„) = 

QK) = (QK)) 

U(u n ) = (%„)) 



(9) 
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V(co n ) = (V(w„)> 

For unpolarized continuum radiation, the harmonic components of different polarization and frequency are independent 
random variables having zero mean and the same variance. Therefore, after taking the ensemble average, the ensemble- 
averaged Stokes parameters Q(u) n ), U(u) n ) and V(u) n ) vanish. We note that the information on the amplitude and relative 
phase of different spherical components of the radiation field at each frequency band is lost after the step of taking the ensemble 
average. That is the fundamental difference between the four parameters (X, Q,M, V) defined above for each realisation and 
the Stokes parameters. 



2.2 Radiation-matter interaction 

The basic transfer equations of maser radiation are derived here for the case involving the transition J=l-0. For higher 
transitions the procedure is similar but algebraic manipulation is considerably more involved. Because we consider a one- 
dimensional maser amplifying the background continuum radiation, the spontaneous emission of masing molecules is ignored. 
This assumption is often used in theoretical studies of astronomical masers. To study self-consistently the effect of spontaneous 
emission, full quantum treatment of the radiation field is necessary and will be the subject of future publications. Goldreich 
et al. (1973), Deguchi & Watson (1990) derived the transfer equation for the radiation field within the framework of the 
rotating-wave approximation. The transfer equation (Goldreich et al. 1973, Deguchi & Watson 1990) can be written as 
follows: 

or in the frequency domain: 

±E v (z,u,)=*™*P v (z,u>) (11) 
az c 

The polarization vector P describes how the masing medium interacts with the radiation field at the frequency of the maser 
line. In this paper we will consider the case where the interaction between the radiation field and the molecule can be described 
within an electric dipole approximation. To calculate the value of the polarization vector P we need to use the density matrix 
p(z,v,t) to describe the medium, which is a collection of masing molecules at a given position z and moving at velocity v. 
The density matrix can be written as: 

/ p++ P+b \ 



p(z,v,t) = 



poo pob 
p— p- h 

\ Pb+ PbO Pb- Pbb / 



(12) 



The Hermitian property of the density matrix implies that pba = p a b where a = ±,0. By writing down the above form of 
the density matrix we have made the assumption that off-diagonal elements p aa / (a^a') are very small, which implies that 
Zeeman splitting is much larger than the stimulated emission rates and thus the mixing of population of different magnetic 
sub- levels does not occur (Goldreich et al. 1973). 

To work in the frequency domain, we first expand the density matrix elements into Fourier series for a given time interval T: 

n — + oc 

n— — oc 

i= +oo 

f>»(z,V,t) = Paa(z,^K^ (t ~ Z/C) (13) 



1— — oo 
i— + 00 



Pbb(z,v,t) = J2 Pbb(*,w,^)e^ ( *- Z/C) 

i— — oc 

Because the diagonal elements of the density matrix p aa and pbb are real, we have the following relations p aa (z, v, cj)=p aa (2, v, —ui) 
and similarly for pbb(z,v,ui). The evolution equation of the density matrix p(z,v,t) can be written in the following compact 
form (Icsevgi & Lamb 1969, Sargent et al. 1974): 



(dt ~*~ V ~diz) Psh = ~( iuJab + r )P ab ~ -^V ab (pbb - p aa ) 
/ d d \ i 

(^+W^-)p aa = A a — T a p aa — -(V ab pb a — V ba pab) (14) 

/ d d \ i 

Iq£ + V ~Q^J P hh = A b - T b pbb - ^(V ba p a b — V ab pb a ) 
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Where w ab = loo + a • gu>B are the frequencies of Zeeman components of the J = 1 — transition, A a (v) and A b (t>) are the 
pumping rates into the upper and lower maser levels. Y a and Y b are the loss rate due to pumping and collisional decoherence 
(Sargent et al. 1974). For the sake of simplicity, we assume here that the loss rates are the same for lower and upper levels of 
the maser transition, v is the velocity of masing molecules. Substituting the Fourier expansion of density matrix elements into 
the above equations (except for the terms involving the interaction matrix V to be written out explicitly later) and collecting 
term by term, we obtain: 

p ah ((J n ,v) = -^-[Vab(/9bb - Paa)(w n )] • ^{iO n ,v) 

Paa(^n,") = j A a (v) 5 n ,0 - ^ [V ab pba ~ V ba p ab ) (w„ )] | • 7+"* (uJ n , v) (15) 
pbb{iO n ,v) = |Ab(u)5n,0 - ^[V b aPab - V ab pba)(^n)]| • T+'K.") 

The 7 functions are the Lorentzian response of the masing molecules to the radiation field and given as follows: 

~/±(iOn,v) = 1/ jr ± i [oj a b - (WO + Wn) • (1 - ^)]| ~ 1/ |r ± i [wo ^ - (wo + W„ - W a b)]j 



7 ±K,w) = V [r t »w„-(i - ^) 



i/[r=F^n] (16) 



7 ^ b K,^) = 1/ r t *w„ • (i - ^) ^ 1/ [r t »w„] 

We note that in arriving at the above equations we have used the similar approximation as in Section II of Menegozzi & Lamb 
(1978) and explained in more detail in Dinh-V-Trung (2009). 

In the frequency domain ui n , the interaction term V • p can be written as the convolution: 

9= + oo 

V-p(w„)= VK- q )-pK) (17) 

q — — oo 

The interaction matrix V between the electric field E and the masing molecules having a dipole moment d can be calculated 
following Deguchi & Watson (1990): 



V= -(E ■ d) = Re|^E M d_ M | =Re j^E p d_ p j ~ I^ Ep d_ p (18) 

In the last step we have used the rotating wave approximation, retaining only the positive frequency part of the electric field 
E. The values of the spherical components d p b of the dipole moment matrix between the upper magnetic sub-level a and lower 
level b of the transition J=l-0 can be easily evaluated as follows (Deguchi & Watson 1990): 

dp= + i = dM=+i(l + cos8)/2 + idM =0 sin8/V2 — i!m=-i(1 — cos6)/2 

dp=_i = -d$=+i(l - cos6)/2 + id$ =o sin0/V2 + d^ = -i(l + cos6)/2 (19) 
We also define here the complex conjugate components of the dipole moment: 

(d£)* = (-l) M d» b M (20) 
The interaction matrix written in the frequency domain has the following form: 

v ab (u;„) = 1 Yl d - b p s p^™) ( 21 ) 
p=± 

Therefore from Eq. 15 we obtain: 

pab(^n,w) = ^2 d - b p £ pK-<l) ' [Pbb(i0q,v) - p a a(^ q ,u)] • 7+ b (cJ„,w) (22) 

p=±, q 

Consequently, Eq. 17 becomes: 

V ab • Pab(^n) = ^ Y (d! b p /)*£p(^ m+n ) ■ £*/(l^ m _ q ) • [p bb *(^ q ) - paa*(^q)] ■ 7- b (^m,w) (23) 

p,p',m,q 

Similar expression for the complex conjugate term: 

V ab * ■ p a bK) = ^ J2 {*?p)**-p>K{<**-*) ■ S p ,(a; m _ q ) ■ [p bb K) - PaaK)] • 7 +Vm, v) (24) 

P,p',m,q 
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The elements of density matrix can then be written as follows: 

Paa(Wn,v) = \ ^ dj b p {&%, ) * Ep (iJm+n) E*, (aJm-q) [pbb H , v) ~ Paa H , ")] (^m , f) 

I mq pp' 

m — n 

mqpp' 
A a (v)*n,o| 7+*H) 

jObb(Wn,v) = -^2^]| ^ d! b p (d! b p /)*gp(g; m+ n)gp/(^m-q) [pbb(^q,«) ~ PLK, f)] T^K, tl) 
a ^mqpp' 

+ X] ( d - b p)*d! b p/£;;(a; m -n)-Bp'(Wm- q ) [Pbb(^q,w) - PaaH, v)] 7+ b (u>m, v) 
mqpp' 

+ A b ( V )5„,o|7+ m K) (25) 
The polarization vector P of the masing medium induced by the radiation field is calculated as: 

oo 

P(z,t)= J dvtr[p(z,v,t)d] (26) 



or written explicitly using the rotating wave approximation: 

1 

2 



/OO 
dv '^2pa.b(z,v,uj„) dp a (27) 
-oo 



Once the equations of the density matrix components are solved in the frequency domain, we can calculate the polarization 

vector of the maser medium using the above expression. 

Using the normalised homogeneous line profile 0± b defined as follows: 

<pf (uj„,v) = — • ^f{u) n ,v) (28) 

7TC 

the transfer equations Eqs. 11 for maser radiation become: 
dE p (z,ui n ) _ 2?r 2 



dz h 

ap'q 



J dv'^(A%)*A%,E p ,(z,u) n - q ) [p^(z,v,w q ) - p hh (z,v,uj q )]<f)f(uj n ,v 



) (29) 



These equations are of great importance in our work because they govern the change of the amplitudes and phases of the 
background radiation through the maser medium once the density matrix components are known. In the next section we will 
use these equations to follow the evolution of the amplitudes and phases of different Fourier components of the radiation field 
through the maser medium 



2.3 First order approximation: radiation transfer equations for Stokes parameters 

In order to derive the standard transfer equations for maser radiation, we need to make two assumptions: the population 
inversion is constant in time (p aa (oj q , v)=0 for w q 7^ 0) and the Fourier components of the radiation field at different frequencies 
are not correlated. In addition, astronomical masers are known to be broadband, i.e the line width due to Doppler broadening 
is much larger than the homogeneous linewidth F and the transition rate due to stimulated emission. Thus we can safely 
ignore the imaginary part of the homogeneous line profile 7J . After straightforward manipulations of equations Eqs. 25 and 
Eqs. 29, and taking the ensemble average, we obtain the following form of the transfer equations for the Stokes parameters 
of the radiation field: 





( /(w) \ 




( A(u) 


B(u>) 





C(u) \ 


( in 
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B(co) 


AH 








QH 




dz 


U{u>) 










ah 





U(u) 






V V(u>) J 




V C(u) 








A(u) ) 




) 



(30) 



The coefficients A(w), B(w) and C(w) are defined as follows: 

A(w) = ^ • 3 B ■ ^ • i [(A P++ + Ap—) (l + cos 2 9) + 2 Ap 00 sin 2 0] 
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B(uj) = ^-35- £ • i[(Ap++ + Ap__) - 2Ap oo ]sin 2 (31) 
C(u>) = ^-35-^ -i[A P++ - Ap__] cos^ 

In the above equations, B is the Einstein coefficient and related to the spontaneous emission probability A by the well-known 
relation B = A ■ c 2 /2hv 3 and 

/ + oo 
dv 4>t h {u),v)[p ++ (v) - pbb(v)] 
-OO 

/ + OO 
dv</£ b (u;,i;)[poo(t>) - Pbb(v)] (32) 
- OO 

/ + OO 
dv<p~ h {uj,v)\p — (u) - Pbb{v)] 
- oo 

where 0* (^i v ) is the real part of the normalised homogeneous line profile <jf^{ui,v). The density matrix representing popu- 
lation in each magnetic sub-level of the masing molecules is determined through the familiar statistical equilibrium equations: 

Fp++{v) = -R++[P++i v ) - Pbb{v)\ + X+(v) 

Tpoo(w) = --Roo[poo(v) - Pbb{v)] + X (v) (33) 
rp— (v) = -R—[p—(v) - pb b (v)] + X-(v) 

Tpbb(w) = R++[p++(v) - Pbb{v)] + Roo[poo(v) - Pbb(w)] + 
R — [p — (v) - pbb(w)] + A b (u) 

where the stimulated emission rate R++, Roo and R are defined as: 

_ 3B v— v . 1 + cos 2 8 r , . ntrt \ sin 2 _ , . . + h . . . 

R ++ = _ ^z^ [ 2 (a;m) + cosei/ ^ m ) + — QK)h, + K,«)A^ 

m 

-Roo = — ^7 ^ sin 2 6[I (cum) — Q(oJm)]'y? h (iOm,v)Auj (34) 

07T * — ' 

m 

3-B v-^r 1 + cos 2 6 r , , . T , . . sin 2 9 _ . . , _k . . . 

fl— = -g^a ^ J (^m) - cos6»l/(o; m ) + — 2— Q(w m )]7r (^m,v)Ao; 

m 

In the above expressions, the functions 7f b (cj,w) are the real part of the homogeneous line profiles ^(ujjV). 
Assuming a small Zeeman splitting and using the fact that the functions ^ b (cj,v) are sharply peaked in comparison to the 
maser linewidth and have the normalisation 7f b (dj m , v)Alo = it, we can easily recover the standard equations derived 
by Goldreich et al. (1973), Deguchi & Watson (1990). Therefore, the standard radiative transfer equations for the Stokes 
parameters follow naturally from our formulation of the astronomical maser. 



3 SIMULATIONS OF MASER AMPLIFICATION 

In our simulation we will choose the parameters appropriate to an astronomical maser. For simplicity, we choose the loss rate 
r = 1 s^ 1 and the normalized pump rates are assumed to have a velocity dispersion a, with the difference in pump rates 
AA(u) = exp(— v 2 /a 2 ). In our simulation we use 600 modes with a frequency resolution of Alu — 0.75 T around the maser 
line covering the range —a to +a in the velocity domain. The actual velocity resolution Aw will depend on the frequency 
of the maser line, as Av = [c/wo]Aw. For a maser line such as the 1612 MHz OH maser, Av is approximately 2 cm s _1 . 
The corresponding value for the velocity dispersion is 200 cm s _1 Although the velocity dispersion of the maser line in our 
simulations is much smaller than in real astronomical masers, the number of frequency modes and the bandwidth are large 
enough, i.e. the bandwidth of 150 s _1 is much greater than the loss rate of 1 s~\ to capture the main features of the broadband 
radiation field produced by the astronomical masers. 

The Zeeman splitting of the upper energy level J = 1 is taken into account explicitly in our simulations. We adopt a value of 
the splitting qujb = 10 Auj or 7.5 s _1 . The adopted splitting is consistent with our starting assumption of the presence of a 
strong B-field in the maser because <?u>b is much larger than the loss rate F = 1 s _1 from the energy levels involved in the 
maser transition. 

Since we deal with only a partially saturated maser, we consider only 10 harmonic components of the density matrix 
paa, bb(<^ q , v) . As shown later, the number of harmonic components is enough to capture the pulsations of the molecular 
population inversion. We generate the background continuum radiation in a simimlar way as Menegozzi & Lamb (1978) using 
random generator RAN2 from Press et al. (1992). The phase of electric field components is random and uniformly distributed 
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over the interval to 2tt. To reduce the enormous amount of computer time required to compute enough realisations of the 
radiation field in order to reach an acceptable level of random noise in the Stokes parameters, we choose to use a constant 
amplitude for all the modes of the background radiation field. Our choice of course allows us to study only the amplification 
of the background radiation field with radom phases. However, the difference in phases between frequency modes of the radi- 
ation field is enough to randomize the Stokes parameters Q and U, which characterize the linear polarization property of the 
background radiation field. Therefore, in our simulations the amplitude of all frequency modes is a constant I±(u> n )ALu — 1 on 
scale of 2hv 3 /c 2 . The intensity of the maser shown in all figures is also of the form 7(w„)Aw, where I(uj n ) is the intensity per 
unit frequency v as defined in the previous section. We use a fourth-order Runge-Kutta method with a fixed step h — 0.02 to 
integrate the transfer equations Eqs. 29. Each realisation of the background radiation field is evolved through the maser by 
solving the equations Eqs. 25 to determine the density matrix, together with integrating the equations Eq. 29 to calculate 
the change of the amplitudes and phases of different spherical components (E- and E + ) of the radiation field. We then record 
the emergent radiation for later analysis. For the sake of simplicity, we will present the results in the form of the parameters 
(X, Q,U, V), as defined in Sec. 2.1 for each realisation of the background continuum radiation. These parameters have the 
advantage of containing the same information on the amplitudes and phases as the spherical components E- and E+, and at 
the same time are directly related to the usual Stokes parameters through a simple ensemble average. 

We carry out our simulations in the partially saturated regime, which is likely relevant to most astronomical masers. The 
choice is also necessary because we consider only a limited number of harmonic components of the molecular population 
inversion. In the saturated regime the fluctuation will be stronger and thus require the consideration of a larger number of 
harmonic components. The amplification of the maser with length L is specified by the unsaturated optical depth t(v = 0) 
at the line center: 



A value of t(v = 0) = 20 is used throughout in our simulations. The spontaneous transition rate between the upper and 
lower maser levels is taken to be 10 -9 s _1 . The Einstein coefficient B is related to the spontaneous transition rate by the 
well-known relation A — [2hv 3 /c 2 ] B. Our choice of unsaturated optical depth corresponds to a partially saturated maser. As 
such, the number of harmonic components of the density matrix used in our simulation is adequate to capture the pulsations 
induced by the radiation field. However, for higher optical depth, a larger number of harmonic components will be necessary. 
We consider two representative cases: in the first case the magnetic field B makes an angle 8 = 30° with respect to the maser 
axis and in the second case the magnetic field is perpendicular (8 = 90°) to the maser axis. The partially saturated masers in 
these two cases allow us to assess the validity of the standard formulation and also test the predictions of the model presented 
in Elitzur (1992, 1993, 1996). 

In Fig. 1 we show the parameters (X, Q, U, V) of one realisation of the background continuum radiation for the case 6 = 30°. 
Because in our simulations the amplitude of both E- and E+ of quasi-monochromatic modes are the same, the parameter V is 
identical to zero. The intensity T is also constant across the bandwidth of the incident radiation. Only the parameters Q and 
U of the radiation field fluctuate strongly due to the random phase of the electric field components E + and E- . However, the 
ensemble average Stokes parameters Q and U are identically zero, as expected from a non-polarized background continuum 
radiation. The amplified radiation at the output of the maser as shown in Fig. 2 clearly displays the line narrowing effect 
and intensity fluctuations induced by the population pulsations. Because the maser is only partially saturated, the harmonic 
components p(z,ui q ,v) of the magnetic sub-level populations are small for ui q ^ as seen in Fig. 3. In Fig. 4 we show the 
ensemble averaged Stokes parameters of the maser radiation calculated using 2400 realisations of the background continuum 
radiation field. These Stokes parameters are consistent with the results of standard formulation for the same case as shown 
in Fig. 5. Similar results for the case 8 = 90° are shown in Fig. 6 and 7. 



4 DISCUSSION 

Our simulation results show that the amplified radiation propagating close to the direction of the magnetic field (8 = 30°) 
possesses both linear and circular polarization. The appearance of circular polarization is a natural consequence of taking into 
account the Zeeman splitting, which in our case is not negligible in comparison to the maser linewidth. Since the stimulated 
emission rate is quite small even at the line centre (R±±/T ~ 0.05, i?oo/r ~ 0.1), the fractional linear polarization Q/I is 
small ~ —1% at the line center, well below the limit of 100% predicted by Goldreich et al. (1973) for fully saturated masers. 
Within the partial saturation regime, which is probably most relevant to astronomical masers, our results are consistent 
with the results obtained from standard formulation (Eqs. 30 and Eqs. 33) and with the calculations of Western & Watson 
(1984), which indicate that linear polarization approaches the 100% limit very slowly for sin 8 ^ 1/3. The fractional linear 
polarization of the maser polarization grows faster in the case 8 = 90 and reaches Q/I ~1.5% at the line center (Figs. 6 & 
7). This result is also consistent with that predicted by Western & Watson (1984). We note that by reducing the multi- level 
model to the idealized two-level case, Gray (2003) also reached similar conclusions. Thus our results agree with previous 




c AA(u) 

' v r 



(35) 
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analytical and numerical simulation works done by Gray (2003) and by Western & Watson (1984). 

The creation of linear polarization when the magnetic field is close to the maser axis (9 = 30° or sin 2 = 1 /4) contradicts 
the model of maser polarization described in Elitzur (1993, 1996), which predicts that propagation of polarized radiation is 
inhibited in the directions close to the magnetic field B. No instability associated with the evolution of fully polarized modes 
of maser radiation is seen during the integration of the radiation transfer equations Eqs. 25 and Eqs. 29. Such instability 
is identified in Elitzur (1993) as the main reason for maser radiation to be linearly polarized at the polarization limit Q/I 
= (3sin 2 9 - 2)/3sin 2 9 for sin 2 9 ^ | even when the maser is unsaturated. Contrary to this assertion, our simulation for 
the case of a partially saturated maser with 9 — 90° shows that the fractional linear polarization of the maser radiation is 
actually small, ~1.5%, well below the limit of 33% given by the above expression. 

We note that in Elitzur (1993, 1996) the equations normally reserved for Stokes parameters and steady-state populations, 
quantities obtained by ensemble averaging over large number of realisations of radiation fields, are used to describe the 
evolution of quasi-monochromatic and fully polarized modes of the radiation field. Consequently, it leads to some difficulty in 
understanding the requirement given in Elitzur (1993, 1996) to perform another averaging step to justify the mathematical 
solutions of the eigenvalue problem for the standard radiation transfer equations. The confusing description of the way to 
handle the random radiation field (Elitzur 1991, 1993) seems to indicate that a self-consistent treatment of the interaction 
between radiation and masing medium has not been achieved. 



5 CONCLUSION 

The explicit incorporation of broadband random radiation field into our treatment of the astronomical maser has proved very 
important to investigate the properties of maser emission. We are able to directly simulate the amplification of the background 
continuum radiation by the maser and study in detail the appearance of linear polarization. Our simulation results show that 
there is no problem with previous numerical studies of the maser polarization in the unsaturated and partially saturated 
regime. Hopefully, the current formulation can be extended and applied to study maser emission under conditions other than 
that considered here. 
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Figure 1. One realisation of the background radiation field. The Stokes parameters (I, Q, U, V) of the quasi- monochromatic modes 
across the bandwith of the radiation field arc shown in each frame. All the parameters (I, Q,U,V) are expressed on a scale of 2hu 3 /c 2 . 
The Stokes parameters I and V have constant values because of the assumption on the constant amplitude of the frequency modes of 
the background continuum radiation. 
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Figure 2. The emergent radiation field amplified by the maser, corresponding to the incident radiation field shown in Fig.l. The magnetic 
field makes an angle 0=30° with respect to the maser axis. The Stokes parameters (X, Q, U, V) of the quasi-monochromatic modes across 
the maser line profile are shown in each frame. All the parameters (I, Q,U, V) are expressed on a scale of 1hv' i j(? . 
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Figure 3. Harmonic components of the magnetic sub-level populations of molecules with zero velocity (v=0) at the end point of the 
maser (z = L). The corresponding emergent radiation field is shown in Figure. 2 
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Figure 4. Stokes parameters of the maser radiation field determined from the ensemble average of the parameters (I, Q, U, V) over 
2400 realisations of the radiation fields. The magnetic field is inclined at an angle = 30° with respect to the maser axis. The Stokes 
parameters are expressed on a scale of t lhv i j<? . 
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Figure 5. Stokes parameters of the maser radiation field predicted using the standard formulation in the case 6 = 30°. The Stokes 
parameters are expressed on a scale of 2hv' i /c 2 . 
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Figure 6. Stokes parameters of the maser radiation field determined from the ensemble average of the parameters (I, Q,U,V) over 
2400 realisations of the radiation field. The magnetic field is inclined at an angle 8 = 90° with respect to the maser axis. The Stokes 
parameters are expressed on a scale of 2hu 3 /c 2 . 
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Figure 7. Stokes parameters of the maser radiation field predicted using the standard formulation in the case 9 = 90°. The Stokes 
parameters are expressed on a scale of 2hv 3 /c 2 . 



